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SIMILARITY  SOLUTIONS  FOR  CONVERGING  SHOCKS 


by 

R.  B.  Lazarus  and  R,  D.  Richtmyer 


ABSTRACT 

This  report  recapitulates  the  known  results  for 
similarity  solutions  for  the  flow  problem  of  a  strong 
converging  shock  in  spherical  or  cylindrical  symmetry 
and  extends  that  work  in  four  ways:  (1)  parameters  of 
the  standard  solutions  are  given  for  a  large  number  of 
values  of  y;  (2)  some  new,  non-analytic  solutions  are 
exhibited  for  relatively  large  values  of  y;  (3)  the 
standard  solutions  are  examined  more  thoroughly  in  the 
limits  and  y-^-l ;  and  (4)  solutions,  existing  only 

in  a  narrow  band  of  values  of  y,  are  given  for  the 
problem  of  two  converging  shocks. 


I .  INTRODUCTION 

As  is  well  known , ^ ^ ^ ^  there  is  a  similarity  solution  for  a 
shock  converging  on  the  origin  in  spherical  or  cylindrical  symmetry, 
when  that  incoming  shock  runs  with  infinite  Mach  number  into  uniform 
material  at  rest  and  when  that  material  obeys  a  gamma  law  equation 
of  state  pe  =  p/(y-l),  with  e  the  internal  energy  per  unit  mass,  p 
the  density,  and  p  the  pressure.  The  solution  includes  reflection 
of  the  shock  at  the  origin,  and  divides  space- time  (r,t)  into  three 
regions,  namely  Region  1  ahead  of  the  incoming  shock.  Region  3  be¬ 
hind  the  reflected  shock,  and  Region  2  between  the  shocks  (see  Fig. 
1-1). 

Previous  authors  have  observed  that  the  solution  is  unique 
(given  gamma  and  the  type  of  symmetry)  if  one  requires  continuity 
of  the  derivatives  of  the  flow  variables  throughout  the  interior  of 
Region  2,  and  the  present  work  includes  calculations  of  those  "stand¬ 
ard"  solutions  for  many  values  of  gamma  (and  both  symmetries),  in¬ 
cluding  the  limiting  cases  y^l  and  y-^".  But  the  present  work  also 
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shows  that  those  solutions  are 
not  unique . 

Even  with  the  requirement 
of  continuous  derivatives,  it  is 
shown  in  Sec.  VI  that,  for  a 
narrow  band  of  values  for  gamma 
near  y  =  2 ,  two  new  solutions 
exist  in  which  Region  2  is  di¬ 
vided  by  a  second  incoming  shock, 
which  overtakes  the  original 
shock  at  the  origin. 

Furthermore,  as  discussed 
in  a  previous  report  in  this 
series,^  for  the  case  of  a  col¬ 
lapsing  cavity,  the  original 
flow  equations  do  not  require 
continuity  of  derivatives.  In 
particular,  there  is  in  Region  2  a 


Fig.  1-1 

r-t  trajectories  of  incoming  and 
reflected  shocks. 


limiting  negative  characteristic 
which  reaches  the  origin  concurrently  with  the  incoming  shock; 
jumps  in  the  derivatives  of  the  flow  quantities  can  be  propagated 
along  that  characteristic.  Since  the  trajectory  of  that  character¬ 
istic  corresponds  to  a  single  value  of  the  similarity  variable,  we 


may  accept  similarity  solutions  with  jumps  of  derivative  at  that 
point.  It  is  shown  in  Sec.  Ill  that  one-parameter  families  of  such 
solutions  exist  for  gammas  greater  than  certain  critical  values  (in 
fact,  the  very  values  which  are  the  threshold  for  the  double-shock 
solution  band) ,  and  that  those  solutions  appear  to  be  quite  inter¬ 


esting. 

Finally,  it  appears  that  these  two  types  of  new  solution  can 
be  combined. 


II.  THE  FLOW  EQUATIONS 

Given  an  inviscid  fluid  without  heat  conduction,  described  by 
its  local  velocity  v(r,t),  its  density  p(r,t),  its  energy  per  unit 
mass  e(r,t),  and  its  pressure  p(r,t),  as  given  by  some  equation  of 
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state  p  =  p(p,e),  the  equations  governing  regions  of  smooth  flow 
are,  in  the  absence  of  any  body  forces, 

Lp  =  0 

Lpv  =  -Vp  (2 .1) 

Lpe  =  -pV*v, 

where  the  operator  L  is  defined  by  Lf  =  f  +  V*vf,  the  subscript  t 
denoting  partial  differentiation  with  respect  to  t.  For  a  poly¬ 
tropic  fluid,  the  equation  of  state  is  p  =  (Y-l)pe,  and  the  entropy 
is  a  function  only  of  the  combination  s  =  pp~^.  Substituting  for  e 
in  the  equations  above,  we  find  that  (v*V  +  9/8t)s  =  0,  so  that  the 
entropy  is  indeed  constant  along  the  trajectories  of  fluid  elements. 

Introducing  the  new  variable  c(r,t)  =  +'^( 9p/ 9 p)  ^  =  +  VyP / P  > 
the  local  sound  speed,  we  can  rewrite  our  equations  as 

+  V* (pv)  =  0 

It 

^t  *  ^  (Y"1)cV*v  =  0. 


For  the  cases  of  cylindrical  symmetry  (v=l)  or  spherical  symmetry 
(v=2) ,  these  can  be  written  (using  u  for  the  radial  fluid  velocity) 

P^  +  (pu)^  +  vpu/r  =  0 

^t  ^^r  ^  ® 

Ct  +  uc^  +  (Y-l)c (u^+vu/r)  =  0. 
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Another  attractive  choice  of  dependent  variables  replaces 
p(r,t)  by  s(r,t)  =  With  the  subs  ti  tut  ion  k  =  2/(y”1)> 

this  choice  yields 

s^  +  us  =0 

t  r  (2.4) 

(u  ±  kc)^  +  (u  ±  c) (u  ±  kc)^  =  +vuc/r  +  c^s^/y (Y“1) s , 

displaying  the  equations  in  characteristic  form. 

Now,  with  a  and  k  free  parameters,  we  try  the  similarity  vari¬ 
able 


y  =  const.  +  log  r  -  alog  t,  (2.5) 

and  the  substitutions 

uCr,t)  =  -  art  ^V(y) 

c(r,t)  =  -  art  ^C(y)  (2.6) 

p(r,t)  =  pQr'^R(y)  or  s(r,t)  =  s^  a^r^t (y)  • 

Using  9f(y)/3r  =  f'/r  and  9f(y)/at  =  -  af’/t,  we  can  substitute 
these  into  Eq.  2.3.  We  find  that  we  get  common  factors  of  p^,  a, 
r,  and  t;  dividing  these  out  and  using  the  more  convenient  X  =  1/a, 
we  derive’ 

R’  +  (k+v+1)RV  +  (RV) '  =  0 

■V(A+’V)  +  V'(l+V)  +  C  [  (k+2)  C+2C '+C '/R]/y  =  0  (2.7) 

2C(X+V)  +  2C'(1+V)  +  (y-1)C((v  +  1)V+V’)  =  0. 

Using  the  variables  u,  c,  and  s,  we  would  get 

(1+V)S'  +  S  [ (2-k(y-1))V+2X]  =  0 
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and,  after  multiplication  by  (1+V+C) , 


[(1+V)2-  c2](V'±kC')  =  vVC[C?(l+V)] 

-  c2[l+C/Cl+V)]  (K+kCA-l))/Y  (2.8) 

-  (l+V+C)  (X+V±C)  (V±kC)  , 

where  of  course  Eq.  2.8  denotes  two  equations,  one  with  all  the 
upper  signs  and  one  with  all  the  lower  signs. 

By  a  bit  of  algebra,  we  can  get  from  the  R,  V,  C  equations  two 
different  expressions  for  1/(1+V),  one  involving  constants  R’/R 
and  V'/Cl+V),  and  the  other  involving  constants  C'/C  and  V'/(l+V). 
Equating  them,  we  can  get  an  expression  whose  derivative  with  re¬ 
spect  to  y  vanishes,  leading  to  a  constant  of  the  motion  and  thus 
reducing  our  system  to  a  system  of  two  equations.  Explicitly,  the 
constant  of  motion  is 

exp  (2y/a)  [R(l+V)  ]*^/R^  ^  =  const.,  with 

q  =  [k(y"1)+2 (l-a)/a]/(K+v+l) .  (2-9) 

With  more  algebra,  we  can  then  put  our  system  into  the  form 


V  =  N^(V,C)/D(V,C) 
C  =  N2(V,C)/D(V,C), 

where 

D(V,C)  =  (1+V)2-C^ 


(2.10) 


(2.11) 


and 


Nj^(V,C)  =  -V(1+V)(X+V)  +  C2  [(v+l)V+2^^2_-<_  ^ 

N2(V,C)  =  -1/2C[V^(2+v(y-1))  +  V( (3- y) X+v (y-l) + 
+Y+1)  +  2X]  +  C3[i  +  <11:^1^21^^]. 


(2.12) 
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Since  the  similarity  variable  y  does  not  appear  explicitly  in 
D,  N  or  N^,  our  system  o£  two  ordinary  first  order  non-linear 
differential  equations  is  autonomous,  and  we  can  write  it  as  a 
single  equation 

dC/dV  =  f (V,C) .  (2.13) 

An  initial  condition  for  this  equation,  however,  is  a  condition  on 
one  branch  of  the  curve  r(t)  =  const. [tl^,  corresponding  to  some 
constant  value  of  y. 

The  free  parameter  k  allows  us  to  handle  the  s  =  constant 
boundary  condition  of  the  cavity  collapse  problem  described  in  Ref. 
5,  by  taking  k  =  -2(X-1)/ (yl) »  and  the  p  =  constant  boundary  con¬ 
dition  of  an  infinitely  strong  shock,  by  taking  k  =  0. 

Since  Eq.  2.13  does  not  contain  y,  it  will  be  convenient  to 
change  similarity  variable  to  x  =  -e'^^* 


III.  THE  CONVERGING  SHOCK;  THE  SINGULARITIES 

To  permit  a  similarity  solution,  any  shock  must  be  at  a  con¬ 
stant  value  of  the  similarity  variable  x  (or  y) ,  so  that  the  physi¬ 
cal  boundary  condition  along  the  shock  trajectory  =  r(t)  can 

be  a  boundary  condition  at  some  x^  for  the  similarity  equations. 

Thus  we  must  have  r  ,  ,  =  constant ‘t*^ .  Note  that,  if  two  or  more 

shock  ,  j-  .  t. 

shocks  exist  within  one  solution,  they  must  all  have  that  form  with 

the  same  value  of  a,  differing  only  through  different  constants  of 

proportionality.  We  are  interested  in  solutions  for  the  range 

0<a<l . 

The  jump  conditions  across  a  shock  become,  in  terms  of  the 
similarity  variables  V,  C,  and  R, 


2C 


R,  =  R(,(1*V„). 


0^ 


(3,1) 
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For  the  initial  shock  converging  into  material  at  rest,  let  us 
take  t  =  0  to  be  the  time  of  shock  collapse,  set  “  ^.(-t)®^, 

and  take  x  =  tCA/r)^  as  our  similarity  variable.  Then  the  shock 
path  is  X  =  -1.  The  solution  for  x  <  -1,  the  initial  and  undis¬ 
turbed  region,  is  simply  u=V=0=c=C,  and  R(x)  =  1  (remember 
that  we  are  now  taking  k  =  0) .  Then  the  jump  conditions  give  us 
the  starting  values 


V(-l)  =  -2/(y+1) 

C(-l)  =  +  V  2yCy-1)  /(Y+1)  C3.2) 

R(-l)  =  (y+1)/(y-1). 

The  solution  must  extend  through  x  =  -0,  which  corresponds  to 
all  of  r  >  0  at  t  =  -0,  and  continue  through  positive  values  of  x 
(and  thus  t)  until  we  get  to  the  reflected  shock.  At  x  =  0,  we 
must  have  V  =  C  =  0,  so  that  u  and  c  may  be  finite  at  finite  r  (see 
Eq.  2.6).  But  then  the  denominator  in  Eq.  2.10,  namely  D  = 
(1+V)^-C^,  will  be  +1,  whereas  it  starts  out  negative  (namely 
-  (Y"1)/ (Y'*’!)  )  ♦  Thus  it  must  pass  through  zero,  but  it  may  do  so 
only  if  the  numerators  and  vanish  simultaneously.  In  other 
words,  our  solution  of  dC/dV  =  must  pass  through  a  singular¬ 

ity  of  the  form  0/0.  It  will  not  do  so  automatically  but  must  be 
made  to  do  so  by  a  suitable  choice  of  the  parameter  a.  Specifi¬ 
cally,  we  will  find  that  a  unique  a(Y)  (for  the  spherical  case,  and 
a  different  unique  a(Y)  for  the  cylindrical  case)  gives  the  ’’stand- 
ard"  smooth  solutions,  but  that,  for  y  large  enough,  other  values 
of  a  give  other  valid  solutions.  To  understand  the  matter,  we  must 
investigate  the  singularities  in  more  detail. 

It  should  first  be  noted  that,  if  we  substitute  =  (1+V)^, 
which  is  to  say  D  =  0,  into  either  or  N2  of  Eq.  2.12,  then  the 
other  N  will  vanish  identically. 
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N, 


Substituting  =  (1+V)^  into  the  first  Eq.  2.12  and  setting 
0  yields  the  cubic 


0 

with 


a 

b 


(1+V) [V2  +  aV  +  b] , 
1  -  rA-l)(Y-2) 


(3.3) 


(3.4) 


The  solution  V  =  -1  is  irrelevant  to  the  converging  shock  problem 
(it  is  the  starting  point  for  the  collapsing  cavity  problem) .  The 
other  two  solutions  are  real  when  the  discriminant 

a2-4b  =  1  -  2^(y+2)  +  (3.5) 

is  positive.  The  discriminant  is  positive  for  \-l  in  the  range 

0<X-1  <  _ yy _ (3.6) 

(Vy+VI)  ^ 

and  it  is  in  that  range  that  we  will  look  for  solutions.  (The  dis- 
criminant  is  again  positive  for  X-1  >  vy/ (VY""'^)  >  this  range  does 

not  seem  to  provide  any  solutions.) 

Note:  For  the  collapsing  cavity  problem,  Eqs.  3.4  through  3.6 

come  out  to  be  the  same  expressions,  but  with  y  replaced  by  Y"1* 
Observe  that  the  two  singularities  are  at  (V,C)  =  (~1,0)  and 
(0,1)  when  X  =  1,  and  move  toward  each  other  as  X  increases.  We 
will  distinguish  the  two  singularities  by  calling  them  "left"  and 
"right"  according  as  we  choose  the  minus  sign  or  the  plus  sign  in 

V  .  =  l/2(-a±Va^-4b) .  (3.7) 

sing 
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It  will  turn  out  that  the  "standard"  solutions  pass  through 
the  left  singularity  for  small  y  and  through  the  right  singularity 
for  large  y.  By  continuity,  then,  there  must  be  critical  values 
for  Y  (one  for  spherical  symmetry  and  one  for  cylindrical) ,  for 
which  the  standard  solution  has  A  at  the  top  of  the  range  given  in 
Eq.  3.6  and  passes  through  the  coalesced  singularity.  It  appears 
that  "non-standard"  solutions  exist  only  for  y  greater  than  these 
critical  values,  which  are 


=  1.9092084,  for  v  =  1  (cylindrical), 
Y^  =  1.8697680,  for  v  =  2  (spherical). 


(3.8) 


For  any  specific  v  and  y»  now,  other  than  one  of  the  critical 
pairs,  let  us  consider  solutions  of  Eq.  2.13  for  some  A  slightly 
displaced  from  the  unique  A(v,y)  which  gives  the  "standard"  solu¬ 
tion.  Consider  the  solution  as  it  approaches  the  singularity 
(which  will,  of  course,  have  been  slightly  displaced  by  the  change 
in  A).  If  the  singularity  is  at  (V^,C^),  say,  we  must  have 

dc  tv-  v.jayav  ^  (c-cpaN^/ac  ^3  ^ 

dV  (V-V  13N,/3V  t  (C-C  )3N,/3c’ 

S'^  1  s  1 


where  the  partial  derivatives  are  evaluated  at  (Vg»C^)  and  are 
simply  algebraic  functions  of  v,  y,  and  A. 

The  general  solution  of  this  equation  is 


[(C-C^)  -  L2(V-V^)]  =  const-[(C-Cp 


L^(V-V^)]  ",  (3.10) 


where,  with 


R 


2 


9N.  9Nt  ^  9N.  9Nt 

-  — £•)  +  4 — £  — i, 

9V^''  ^9V-  9C^’ 

s  s  s  s 


(3.11) 
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we  have 


(3.12) 


2Li^28Ni/8Cs  =  aN2/9C^-3N^/9V^±R 

and 

2E^  "  3N2/9Cg+9N^/9Vg±R.  (3.13) 

For  our  case,  it  appears  that  R  is  always  real  and  non- zero,  and 
that  and  L2  have  opposite  signs,  in  a  neighborhood  of  the  "stand¬ 
ard"  A(v,y)*  The  E’s  and  L's  are  of  course  algebraic  functions  of 

V ,  Y  >  A  . 

If  E^  and  E^  have  opposite  signs,  the  only  solutions  through 
the  singularity  are  (locally)  the  special  solutions 

C-Cs  -  Lj_2(V-V3). 

For  Y  <  Y  ,  the  standard  solution  is  of  this  type.  For  one  particular 
value  of  A,  the  solution  passes  through  the  left  singularity  with 
the  slope  corresponding  to  the  negative  L,  and,  for  that  A,  the  left 
singularity  does  indeed  have  E's  of  opposite  sign.  For  neighboring 
values  of  A,  the  solution  will  not  pass  through  either  the  left  or 
the  right  singularity. 

Nor  does  it  seem  likely  that,  for  y  <  y^,  there  are  other  solu¬ 
tions  for  substantially  different  values  of  A.  For  larger  values, 
the  E’s  continue  to  have  opposite  signs.  For  substantially  smaller 
values  of  A,  the  E's  do  have  the  same  sign,  but  the  left  singularity 
moves  further  to  the  left,  the  positive  L  is  less  than  one,  and  the 
solution  hits  the  forbidden  line  C  =  1+V  before  it  can  be  attracted 

to  the  singularity  (see  Fig.  3-1). 

For  Y  >  Y  ,  where  the  standard  solution  goes  through  the  right 
'  c  ^ 

singularity,  we  have  the  case  where  the  E's  have  the  same  sign.  In 
such  a  case,  all  solutions  which  come  sufficiently  close  to  (^g>Cg) 
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Attractive  singularity  blocked  Attractive  singularity  open  above, 
above  by  the  line  D=0. 


pass  through  the  singularity,  and  they  do  so,  in  general,  asymptot¬ 
ically  like 


C-Cs  =  L.(V-Vg),  (3.15) 

where  (i  =  1  or  2)  is  the  E  o£  lesser  magnitude.  It  turns  out 
that  that  is  the  positive  L,  and  that  the  (unique)  ’’standard" 
solution  is  precisely  the  special  solution  which  goes  through  with 
negative  slope  (i.e.,  with  the  other  L) . 

For  the  entire  range  Y  >  Y^,  the  positive  L  is  greater  than 
one.  To  reach  the  singularity  without  first  crossing  C  =  1+V, 
therefore,  the  solution  must  come  in  from  above  (see  Fig.  3-2). 

Since  the  main  effect  of  changing  X  is  to  move  the  singularity 
(i.e.,  the  solution  curve  does  not  change  much  until  we  approach  the 
singularity),  this  means  that  the  right  singularity,  with  which  we 
are  here  concerned,  must  be  moved  left .  Thus  only  values  of  X 


greater  than  the  standard  X(v,y)  will  work.  The  foregoing  analysis 
is  only  valid  in  a  neighborhood  of  the  singularity.  A  complete 
analysis  will  be  published  elsewhere. 
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Fig.  3-3 

Pressure  profiles  for  y  =  3, 
spherical  symmetry. 

Figs.  3-3  through  3-5  show 
the  pressure,  density,  and  veloc¬ 
ity  profiles  at  a  time  when  the 
incoming  shock  is  at  r  =  1,  for 
the  case  y  =  3,  v  =  2.  For  the 
non-standard  solutions,  the  cor¬ 
ners  are  on  the  limiting  charac¬ 
teristic  and  are  such  as  to  sat¬ 
isfy  the  flow  equations  from  the 
left  and  from  the  right.  Note: 

In  these  solutions,  the  curves  in 
the  V-C  plane  were  allowed  to 
leave  the  singularity  in  the 
"standard"  direction.  This  is 
not  necessary  (see  p.  10  of  Ref. 
5) ,  but  other  solutions  have  not 
yet  been  studied. 

IV.  THE  REFLECTED  SHOCK 

The  initial  shock, 
collapse  is  at  X  =  0; 
continuation  to  positive 


Fig.  3-4 

Compression  profiles  for  y  =  3, 
spherical  symmetry. 


Fig.  3-5 

Fluid  velocity  profiles  for  y  = 

3,  spherical  symmetry.  The  down¬ 
ward  sloping  portions  are  dis¬ 
tinct  but  too  close  to  plot 
separately . 


which  is  our  starting  point,  is  at  x  =  -1; 
and  continuation  to  positive  times  is  simply 
X.  As  one  might  expect  on  physical  grounds 
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it  will  not  be  possible  to  continue  the  same  solution  to  x  =  +  o®. 

One  expects  this  because  large  positive  values  of  x  correspond  to 
small  values  of  r  at  large  positive  values  of  t,  and  this  region  of 
the  flow  should  be  behind  a  reflected  shock.  As  mentioned  above, 
the  trajectory  of  that  reflected  shock  will  have  to  lie  on  x  =  con¬ 
stant  =  3,  say,  for  some  3  >  0. 

If  we  can  find  the  separate  similarity  solution  for  the  region 

A  /N  A 

behind  this  reflected  shock,  say  V,  C,  and  R,  then  the  two  solutions 
will  have  to  satisfy  the  jump  conditions  at  x  =  3*  We  will  need  to 
satisfy 


(4.1) 


6^(3)  =  c^(3) 


+  X^[(l+V(3))^  -  Cl+V(3))^] 


Note  that  the  constant  of  motion  will  be  a  different  constant  on 
the  two  sides  of  the  reflected  shock,  just  as  it  is  a  different  con¬ 
stant  on  the  two  sides  of  the  initial  shock. 

This  separate  solution  is  needed  for  3  ^  x  <  «>,  and  the  only 
thing  we  have  to  serve  as  a  boundary  condition  is  the  following. 

We  want  u(r=0,t>0)  to  be  zero,  by  isotropy,  and  we  want  c(r=0,t>0) 
to  be  finite.  In  fact,  we  expect  u  to  be  proportional  to  r,  for 
small  r,  so  we  expect  V  to  be  constant  and  C  to  become  infinite  as 
x-»-«>. 

The  standard  trick  is  to  take  a  new  variable  w  =  with  a  a 

positive  number  to  be  determined,  and  to  try 


VCkw)  =  Vp  +  V^kw  +  V^Ckw)"'  +  ... 
C(kw)  =-(kw)"^  +  +  C2kw  +  ..., 


(4.2) 


where  k  is  a  free  parameter  (our  differential  equations  are  homo- 

A  A  A  A 

geneous  in  w,  so  that  V(kw),  C (kw)  are  solutions  whenever  V(w),  C(w) 
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are).  Matching  powers  of  kw,  we  find  that,  if  we  take 


Aa 


1  +  Cv+1)  CX-1) 

^  XWi)Y-^CX--i7' 


Vq  =  -2(X-1)/y(v+1) , 


(4.3) 


then  we  get  a  solution. 

If  we  now  think  of  the  jump  conditions  as  an  operator  which  can 
be  applied  to  our  original,  incoming  shock  solution,  for  arbitrary 
positive  X,  then  we  have  "target"  functions  V  (x)  and  (x) ,  to  be 
matched  by  V(kw)  and  C(kw).  The  value  of  x  at  which  that  match 
occurs  is,  of  course,  just  3.  If  the  value  of  kw  at  which  the  match 
occurs,  is,  say,  z,  then  we  can  determine  k  by  setting  z  =  k3“‘^,  and 
we  have  the  complete  solution. 

V.  THE  LIMITS  y^“  AND  y^I • 

For  Y"^"” >  need  only  switch  to  V  =  Y^>  and  then  we  can  go  to 

2 

the  limit  explicitly.  The  denominator  D  becomes  simply  1-C  ;  the 
numerator  for  V’  becomes 

=  -  XV  +  c2[(v+l)V  +  2(X-1)],  (5.1) 

and  the  numerator  for  C’  becomes 

=  C*[C2  -  X  - (v+l-X)V/2] .  (5.2) 

These  somewhat  reduced  equations  can  be  integrated  numerically  by 
the  methods  described  below  for  general  y. 

For  Y-^lj  the  situation  is  slightly  more  complicated,  because 
the  singularity  approaches  the  starting  point  (V,C)  =  (-1,0).  If  we 
define  =  y-1,  then,  to  lowest  order,  our  starting  point  is 

(Vo,Co)  =  (-l+e^/2,e/V2) .  The  starting  value  for  D  is  then  -e^/2. 
Now  if  we  tentatively  assume  that  X-1  will  turn  out  to  be  of  order 
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e,  we  find  that  the  leading  terms  in  are 

-(v+l)c2  +  (1+V)2  +  (X-1)C1+V),  (5.3) 

and  the  leading  terms  in  are 

%  C[C^-(1+V)^]  (X+V)/C1+V),  (5.4) 

considering  that  we  must  integrate  from  1+V  =  until  1+V  =  C  at 
the  singularity. 

If  we  integrate  dC/dV  =  holding  C  fixed  on  the  right  hand 

side,  we  find,  consistently,  that  C  changes  only  by  a  factor  1  - 
order  (eloge) ,  and  we  find,  again  consistently,  that  we  must  have 


X-1  =  V  V(y-1)/2.  (5.5) 

This  is  confirmed  numerically,  as  well  as  the  additional  result 
that  the  Mach  number  of  the  reflected  shock  is  '^2/  (yl)  >  independent 
of  V  (see  Table  5) . 


TABLE  5 


BEHAVIOR 

OF  THE  SIMILARITY 

SOLUTIONS 

AS  GAMMA 

APPROACHES 

UNITY 

(l-a)2 

Y-1 

(y- 

-1)M^ 

(y- 

i)e 

Y-1 

v=2 

v=l 

v=2 

V=1 

v  =  2 

v=l 

0.1 

0.4163 

0.1317 

1.262 

1.685 

1.615 

1.354 

0.01 

0.9630 

0.2755 

2.117 

1.986 

0.001 

1.4864 

0.3953 

2.202 

2.153 

1.518 

1.252 

10'4 

1.7909 

0.4587 

2.102 

2.073 

1.215 

1.109 

lO'S 

1.9216 

0.4849 

2.042 

2.030 

1.084 

1.044 

10"^ 

1.9641 

0.4933 

2.022 

2.014 

1.033 

1.017 

0 

2 

1/2 

2?  • 

2? 

1? 

1? 
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VI.  MULTIPLE  CONVERGING  SHOCKS 

I£  there  is  a  similarity  solution  corresponding  to  more  than 

one  incoming  shock,  then  the  shocks  must  have  the  trajectories 


r . 
1 


A^C”t) 


a 

9 


(6.1) 


with  A  =  <  . . .  I£  be  the  value  o£  the  similarity  vari¬ 

able  on  which  the  Uh  shock  exists,  then  we  must  have  =  -(A/A^)2. 

Consider  i=2. 

With  D(V,C)  =  (1+V)^  -  C^,  the  jump  conditions  o£  Eq.  3.1  imply 


so  that  D  must  change  sign.  Furthermore,  since  Region  2  is  behind 
the  shock,  we  have  >  0,  and  thus  the  third  jump  condition 

implies  (I+V2) ^  <  (1+V^)^.  But  then  the  second  jump  condition  im¬ 
plies  <  C2^,  giving  >  D2.  Since  and  D2  are  o£  opposite 

sign,  we  have  D2  <  0,  >  0. 

Note:  A  "shock"  existing  right  the  singularity  D  =  0  has 

Mach  number  unity  and  is  not  a  shock  at  all. 

Since  our  solution  behind  the  initial  shock  starts  out  with  D 
negative,  we  see  that  x^  must  be  greater  than  the  value  £or  which 
the  region  1  solution  crosses  the  singularity.  Thus  we  must  have 
the  same  value  o£  a  (=1/A)  as  we  have  £or  the  single  shock  case, 
since  a  is  determined  precisely  by  the  necessity  o£  passing  through 
the  singularity. 

Rewriting  the  jump  condition  on  D  in  the  £orm 


“2  = 


D, 


<  D 


1’ 


-^1 

<^1 

(6.3) 

[  Y+1  \  \ 

)y 

that 

D  >  0  implies 

(C^/(l+V^))2 

<  1 ,  we  see  also  that 

with 

jL 

the  inequality 

stronger  £or 

smaller  values  o£  y. 
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This  means  that  the  vector  in 


(  Vq.Co) 


the  V-C  plane  connecting  (V^,C^) 
to  (¥^,€2)  has  negative  slope 
between  -1  and  0. 

When  the  matter  is  investi¬ 
gated  numerically,  it  turns  out 
that  the  locus  o£  points 

as  X2  ranges  toward  zero  from 
,  the  value  of  x  corresponding  to 

the  singularity,  is  an  arc  con¬ 
necting  the  singularity  to  the 
starting  point  (¥^,0^)  and  lying 
always  below  and  to  the  left  of  Fig.  6-1 

the  original  solution  curve  (see  dashed  line  is  the  locus  of 

points  accessible  from  s,  by  the 
Fig.  6-1).  When  an  attempt  is  jump  conditions.  ^ 

made,  however,  to  continue  the 

solution  from  any  of  those  points  (V  ,C  ),  it  develops  that  the 
solution  moves  almost  parallel  to  the  original  solution  curve. 

Hence,  the  continued  solution  cannot  pass  again  through  the  same 
singularity . 

This  immediately  suggests  that  when  y  is  greater  than  the 
critical  value  of  Eq.  3.8,  so  that  the  primary  solution  goes  through 
the  right  hand  (upper)  singularity,  a  point  CV2,C2)  can  be  found  so 
that  the  continued  solution  will  pass  through  the  left  hand  (lower) 
singularity.  This  turns  out  indeed  to  be  the  case  when  y  is  greater 
than  y^  by  an  amount  small  enough  that  the  width  of  the  locus 
(measured  parallel  to  the  45°  line  C  =  1+V)  is  not  less  than  the 
spacing  between  the  two  singularities.  In  fact,  there  will  be  two 
double  shock  solutions,  for  a  band  of  y  values,  corresponding  to 
relatively  weak  and  relatively  strong  second  shocks,  with  the  two 

I 

fi  solutions  coalescing  at  the  top  of  the  band  and  then  ceasing  to 

t  exist  as  y  leaves  the  band. 

I  For  y*s  above  this  band,  non-standard  solutions  may  exist  with 

y's  sufficiently  close  to  the  upper  bound  of  Eq.  3.6,  which  is  to 
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say  with  the  two  singularities  sufficiently  close,  to  permit  fur¬ 
ther  solutions  with  two  incoming  shocks.  A  complete  analysis  will 
be  published  elsewhere. 

The  entire  situation  can  be  grasped  most  simply  as  follows. 
Pick  values  for  v  and  y,  with  y  >  y^ .  This  determines  a  starting 
point  in  the  V-C  plane,  and  a  one  parameter  family  of  (incomplete) 
solutions  labeled  by  X.  Pick  a  value  for  X  which  lets  the  solution 
pass  through  the  right  singularity  and  continue  to  the  origin;  call 
the  corresponding  solution  curve  Now  that  we  have  X,  we  can 

locate  the  unused  left  singularity  and  construct  the  solution  curve 
(call  it  S^)  which  passes  through  it  in  the  standard  direction. 
Lastly,  we  draw  in  the  {N 2,(^2)  locus  corresponding  to  potential 
second  shocks.  Then  we  have  zero,  one,  or  two  double  shock  solu¬ 
tions  according  as  that  locus  cuts  in  zero,  one,  or  two  points, 
because  we  have  a  physical  method  of  jumping  from  solution  curve 
to  solution  curve  S^.  Finally,  if  the  left  singularity  should  have 
eigenvalues  of  the  same  sign,  then  there  would  by  a  family  of  S^’s, 
all  valid. 

Typical  solutions  are  shown  in  Table  6-1. 

TABLE  6-1 

MACH  NUMBERS  FOR  WEAK  AND  STRONG  SECOND  SHOCKS 


For 

V  =  1  (Yp  = 

1.9092084) 

Y 

Ml 

M2 

1.91 

1.000566 

1233.532 

1.95 

1.078474 

22.34995 

2.00 

1.248826 

8.846871 

2.05 

1.548093 

4.701998 

2.10 

none 

For 

V  =  2  (Y^  = 

1.8697680) 

1.87 

1.000170 

3508.718 

1.90 

1.070279 

24.19837 

2.00 

1.848820 

3.496177 

2.009 

2.220365 

2.701936 

2.01 

none 
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VII.  CONDITIONS  BEHIND  THE  REFLECTED  SHOCK 

The  position  o£  the  reflected  shock,  as  a  function  of  time,  is 

given  by 


r  (t)  =  A3“°‘t°‘,  (7.1) 

X  •  o  • 

where  3  is  the  value  of  the  similarity  variable  x  corresponding  to 
the  shock  trajectory,  as  mentioned  above,  and  A  is  the  constant 
appearing  in  the  trajectory  of  the  initial  shock. 

For  the  region  behind  the  reflected  shock  (inside  it,  geomet¬ 
rically  speaking) ,  it  is  of  interest  to  consider  the  time  depend¬ 
ence  of  the  volume  integrals  of  mass,  internal  energy,  and  kinetic 
energy,  and  of  the  "mean  free  path"  integral  of  pdr.  By  appealing 
to  the  original  substitutions  (Eq.  2.6)  for  u,  c,  and  p,  and  by 
substituting  for  r  the  appropriate  expression  in  terms  of  ;t  and  the 
similarity  variable  w  (which  runs  from  zero  to  3  '^)  ,  we  find  the 
following,  for  given  v  and  y,  taking  p^  =  1. 

The  total  mass  is  simply  proportional  to  the  total  volume,  with 
no  other  time  dependence,  and  the  integral  of  pdr  is  simply  propor¬ 
tional  to  r  .  (The  volume,  of  course,  is  going  like 
The  total  internal  energy  and  the  total  kinetic  energy  are  separate¬ 
ly  proportional  to  the  volume  times  the  factor  t’^^^  As  re¬ 

quired  for  physicality,  a  is  always  less  than  unity,  so  that  the 
average  values  of  internal  and  kinetic  energies  per  unit  volume 
decrease  with  time.  (These  results  also  imply  that  the  energy  den¬ 
sities  behind  the  reflected  shock  are  instantaneously  infinite  at 
collapse  time.  This  is  in  accord  with  the  fact  that  C(x)/x  and 
V(x)/x  remain  finite  at  x  =  -0,  so  that  the  fluid  velocity  u(r,t) 
and  the  sound  speed  c(r,t)  behind  the  initial  shock  become  infinite 
like  at  collapse.) 

The  various  constants  of  proportionality  are  given,  as  func¬ 
tions  of  V  and  y,  in  Tables  7-1  and  7-2.  Ij^  and  I2  are,  respec¬ 
tively,  the  internal  and  kinetic  energies  per  unit  volume,  times 
A-232at2 (1-a) ^  I^  is  the  mass  per  unit  volume;  I^  is  the  mass  per 
unit  area  divided  by  r 

r .  s . 
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TABLE  7-1 

VARIOUS  PARAMETERS  OF  THE  STANDARD  SIMILARITY  SOLUTION,  AS 
FUNCTION  OF  GAMMA;  SPHERICAL  CONVERGENCE. 


It 

2  (spherical) 

ct  6 

Mach  # 

P2(3) 

^1 

^2 

^3 

^4 

1.1 

.79596980 

16,1541 

3.5530 

333.6 

1.95+6 

243.33 

4.8+4 

3.6+4 

1.2 

.75714179 

6.43123 

2.3502 

49.05 

1.37+4 

11,0163 

1710. 

1253. 

1.3 

.73377673 

3.81021 

1.9384 

19.03 

1031.6 

2.2220 

320.4 

231.1 

1.4 

.71717450 

2.68849 

1.7356 

10.72 

197.5 

.8029 

114.1 

81.83 

1.5 

.70442807 

2.08773 

1.61553 

7.226 

60.76 

,3894 

56.12 

40.21 

1.6 

.69418951 

1.72065 

1.53621 

5.479 

24.73 

,2243 

33.30 

23.91 

5/3 

.68837682 

1.54790 

1.4967 

4.719 

15.09 

.1655 

25.20 

18.14 

1.7 

.68571652 

1.47617 

1.47982 

4.419 

12.07 

.1443 

22.27 

16.05 

oo 

.67855370 

1.30320 

1.43758 

3.729 

6.685 

.1002 

16.17 

11.71 

1.9 

.67240014 

1.17523 

1 .40470 

3.251 

4.058 

.07355 

12.45 

9.061 

2.0 

.66704607 

1.07725 

1.37833 

2,902 

2.641 

.05628 

10.01 

7.327 

2.2 

.65816533 

.938224 

1.33851 

2.432 

1.3015 

.03609 

7.102 

5,262 

2.4 

.65108461 

.845319 

1.30973 

2.134 

.7399 

.02521 

5.489 

4.115 

2.6 

.64530018 

.779560 

1.28784 

1.930 

.4645 

.01870 

4.493 

3,406 

2.8 

.64048378 

.731006 

1.27056 

1.783 

.3138 

.01453 

3,835 

2.936 

3.0 

.63641060 

.693969 

1.25649 

1.672 

.2232 

.00159 

3.360 

2.600 

3.2 

.63292118 

.664976 

1.24482 

1.587 

.1663 

.009571 

3.024 

2.358 

3.4 

.62989873 

.641817 

1.24487 

1.518 

.1279 

.008029 

2.762 

2.171 

3.6 

.62725578 

.622959 

1.22641 

1.463 

.1008 

.006819 

2.550 

2.023 

3.8 

.62492541 

.607422 

1.21895 

1,417 

.08034 

.005740 

2.354 

1.896 

4.0 

.62285554 

.594419 

1.21246 

1.3785 

.06746 

.005209 

2.263 

1.815 

4.5 

.61857036 

.569946 

1.19897 

1.3056 

.04302 

.003673 

1.960 

1.6Z2 

5.0 

.61522398 

.552999 

1.18871 

1.254 

.03095 

.002974 

1.832 

1.518 

5.5 

.61253956 

.540843 

1.18033 

1.217 

.02222 

.002258 

1.665 

1,419 

6.0 

.61033915 

. 531820 

1.17342 

1.188 

.01773 

.001978 

1.628 

1.374 

6.5 

.60850311 

.524919 

1.16776 

1.165 

.01348 

.001544 

1.503 

1.308 

7.0 

.60694820 

,519578 

1.16281 

1.147 

.01116 

.001359 

1,473 

1.279 

8.0 

.60445829 

.511963 

1.15480 

1.120 

.007824 

.001033 

1.398 

1.227 

10 

.60104880 

.503479 

1.14368 

1.087 

,004273 

.000615 

1,263 

1.149 

50 

.59073010 

.494072 

1.10718 

1.0118 

.3173* 

.06482* 

1.043 

1.018 

100 

.58950281 

.494977 

1.10233 

1.0055 

.2859* 

.05965* 

.9883 

.9975 

00 

.58828929 

.496368 

1,09753 

1 

0.2880* 

.6456* 

*:  (y+1)^I 
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TABLE  7-2 


VARIOUS  PARAMETERS  OF  THE  STANDARD  SIMILARITY  SOLUTION, 
AS  FUNCTIONS  OF  GAMMA;  CYLINDRICAL  CONVERGENCE. 


V  =  1  (cylindrical) 


Y 

a 

6 

Mach  # 

P2(e) 

^1 

^2 

I3 

I4 

1.1 

.88524806 

13.5364 

4.10488 

29,55 

3.28+5 

15.614 

5236. 

4681. 

1.2 

.86116303 

6.09996 

2.78911 

11.15 

6889.0 

1.9614 

506.1 

445.0 

1.4 

.83532320 

2.81561 

2.02295 

4.796 

203.46 

.2822 

65.48 

57.11 

5/3 

.81562490 

1.69479 

1.69965 

2.928 

21.144 

.07873 

19.31 

16.88 

1.8 

.80859994 

1.44082 

1.61796 

2.527 

10.139 

.05152 

13.32 

11.69 

1.9 

.80409908 

1.30515 

1.57247 

2,316 

6.4327 

.03949 

10.68 

9.390 

2.0 

.80011235 

1.19963 

1.53602 

2.154 

4.3370 

.03129 

8.870 

7.817 

2.4 

.78776900 

.941829 

1.44206 

1,763 

1.3299 

.01528 

5.270 

4.694 

3.0 

.77566662 

.763158 

1.37121 

1.496 

.4265 

.007407 

3.398 

3.069 

3.4 

.77000368 

.697702 

1.34349 

1,399 

.2485 

.005175 

2.828 

2.575 

4.0 

.76363465 

.634863 

1.31564 

1.306 

.1329 

.003377 

2.341 

2.150 

5.0 

.75640105 

.575038 

1.28751 

1.219 

.06190 

.001948 

1.917 

1.783 

6.0 

.75156168 

.540788 

1.27050 

1.169 

.03523 

.001273 

1.693 

1.589 

10 

.74182593 

.483613 

1 .23980 

1.0867 

.008622 

.0003962 

1.325 

1.281 

50 

.73002154 

.431537 

1.20756 

1,0142 

.6287* 

.04112* 

1.0736 

1.0554 

100 

.72853594 

.426147 

1.20374 

1,0069 

.5831* 

.03909* 

1.0358 

1.0269 

00 

.727048052 

.421009 

1.199865 

1 

.5431* 

.03745* 

*:  Cy+I)^! 


Another  integral  o£  possible  interest  is  the  integral,  from 
time  zero  to  time  t,  of  the  volume  integral  of  a  power  of  the  tem¬ 
perature  (or  pressure  or  internal  energy;  we  are  dealing  with  a 
polytropic  fluid)  times  some  function  of  the  density.  One  might 
imagine  such  an  integral  measuring  the  total  amount,  taking  place 
up  to  time  t,  of  some  reactive  process  having  such  a  dependence  on 
density  and  temperature  (assuming,  of  course,  that  the  energetics 
of  the  process  do  not  break  the  similarity  solution).  For  the  nth 
power  of  the  temperature,  we  find  the  following  rather  curious  re¬ 
sult.  If  2n  <  (va+a+1)/ (1-a) ,  then  the  integral  is  entirely  regu¬ 
lar  and  goes  like  .  But  for  any  larger  value  of 
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n,  the  integral  would  diverge  unless  other  effects  (such  as  deple¬ 
tion  of  the  reactants)  were  taken  into  account.  For  v  =  2,  the 
critical  values  for  nCy)  for  example,  n(1.4)  =  5.57,  n(3) 

=  4.00,  n(“)  =  3.36. 

VIII.  THE  NUMERICAL  INTEGRATION 

All  calculations  were  done  on  the  Maniac  II  computer,  using 
the  Madcap  V  system.  All  constants  and  variables  entering  into  the 
integration  of  the  differential  equations  were  carried  with  at 
least  16  decimal  digits.  Explicit  fifth  order  Runge-Kutta  was 
used,  with  step  size  controls  as  discussed  below. 

For  the  a  search,  and  in  fact  for  all  the  region  behind  the 
incoming  shock,  the  independent  variable  used  was  x  =  t(A/r)  ,  and 
the  dependent  variables  were  v(x)  =  -V(x)/x  and  c(x)  =  C(x)/x,  The 
minus  sign  is  historical  accident;  the  division  by  x  is  to  give 
nice  behavior  at  the  star  point  singularity  V  =  C  =  0.  The  initial 
value  for  x  is  -1,  and  integration  must  be  continued  past  the  un¬ 
known  value  X  =  3.  An  efficient  method  of  coping  with  this  diffi¬ 
culty  is  described  below. 

The  equations  were  used  in  the  form  dv/dx  =  N^^/D,  dc/dx  =  N2/D, 
where  now 


D  =  (l-vx)2  _  (cx)2^  (8.1) 

Ni  =  p^[v^Cl-vx)  +  p^c^]  -  p^vc^x, 

2  P2 

N2  =  C[v(p3-P20vx)  +  p^c  x(l  - 


and  the  constants  are 
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=  1-a 

P2  =  2/y 

Pj  =  (v+l)a  -  1  (8.2) 

P5  =  1/2  [CY+l)(l-a)  -  av(Y-l)] 

P2Q  =  l-a-av(Y-l)/2. 

For  the  search  for  the  "standard"  a(v,Y)>  we  exploit  the  facts 
that  the  correct  solution  goes  quite  smoothly  (in  the  V-C  plane) 

from  its  starting  point  through  the  singularity,  that  the  positions 

of  the  singularities  are  quite  sensitive  to  the  value  of  a,  and 
that  the  solution  curve  for  a  wrong  value  of  a  does  not  differ  much 

from  the  correct  solution  curve  all  the  way  up  to  a  point  where  we 

can  determine  that  we  do  indeed  have  a  wrong  value.  Accordingly, 
an  efficient  iterative  algorithm  is  to  choose  the  next  guess  for  a 
so  as  to  move  the  relevant  singularity  on  to  the  line  connecting 
the  initial  (v,c)  point  to  the  last  (v,c)  point  reached  before  the 
aforesaid  determination.  In  practice,  this  determination  was  made 
if  dV/dx  changed  sign  or  if  [dv/dx]  became  larger  than  three  times 
its  initial  value.  (When  calculating  the  non-standard  solutions 
discussed  on  pp .  10-12,  the  "determination"  is  simply  suppressed.) 

It  is  important  to  note  that  all  finite  numerical  representa¬ 
tions  of  a  will  be  determined  to  be  wrong  if  we  approach  the 
singularity  with  a  sufficiently  small  step  size.  Conversely, 
almost  any  value  for  a  will  get  us  through  the  singularity  without 
such  determination  if  we  approach  with  a  sufficiently  large  step 
size.  Accordingly,  the  step  size  was  automatically  reduced  to  a 
prescribed  as  we  approached  the  singularity,  and  no  further, 

with  h  .  chosen  to  give  the  desired  accuracy.  The  bulk  of  the  a 

search  work  was  done  with  h  .  =  ^  lO"^. 

mill 

The  code  was  run  in  the  a  search  mode  for  all  desired  values 
of  (v,y),  without  continuing  the  solution  past  the  singularity. 
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Then  it  was  rerun  with  the  correct  a  and  with  larger  (usually 

2"^^  ^  6.4x10“^)  to  get  the  complete  solution.  In  this  mode,  it 
was  almost  always  true  that  the  solution  would  step  smoothly 
through  the  singularity  in  a  single  step,  but  this  is  a  matter 

o£  luck.  I£,  as  happened  occasionally,  the  code  determined  that  a 
step  (or  a  partial  step  within  the  Runge-Kutta)  might  accidentally 
land  too  close  to  the  singularity,  then  it  took  a  "jump”  step  o£ 

8h  .  and  printed  a  noti£ication.  (Note  that,  as  discussed  else¬ 
where,  we  are  passing  through  the  singularity  in  an  eigendirection , 
without  change  o£  slope.) 

In  continuing  the  solution  through  x  =  0,  the  step  size  was 

again  reduced  to  h  .  ,  in  order  to  permit  printing  out  accurate 
^  min’ 

values  o£  V(x)/x  and  C(x)/x  at  x  =  0. 

Two  problems  arise,  now,  in  connection  with  continuing  this 
phase  o£  the  solution  up  to  x  =  3.  We  must  be  sure  to  go  £ar 
enough,  but  we  do  not  want  to  waste  time  going  too  £ar,  and  we  need 
£inely  spaced  tables  o£  the  "target"  £unctions  de£ined  on  page  14, 
but  only  in  the  neighborhood  o£  (the  unknown)  3*  The  two  problems 
are  solved  as  £ollows. 

The  code  is  given  a  lower  bound  £or  3,  call  it  i£  no 

better  in£ormation  is  available,  then  zero  is  the  lower  bound  used 
by  de£ault,  but  o£  course  we  can  do  much  better  than  that  once  we 
have  sketched  out  3(v,y)  by  running  a  £ew  cases.  The  code  then 
saves  the  solution  £or  some  value  o£  x  near  B^^^  and  pushes  ahead 
using  large  steps  and  saving  a  coarse  table  o£  the  target  £unctions. 
It  pushes  ahead  until  it  approaches  the  singularity  C  =  - (1+V) , 
which  must  always  lie  beyond  3>  and  then  £inds  the  (approximate) 
re£lected  shock  solution  (see  below)  and  an  approximate  value  £or 
3.  Then  it  picks  up  the  saved  solution  £rom  near  B^^^  and  moves 
ahead  with  £ine  steps  until  x  is  sa£ely  beyond  the  approximate  3, 
and,  £inally,  gets  an  accurate  solution  £or  the  re£lected  shock. 

For  the  region  behind  the  re£lected  shock,  the  independent 
variable  used  was  t  =  kx"'^,  where  k  is  a  £ree  parameter  that 
cancels  out  o£  the  di££erential  equations  and  is  used  as  described 
below,  and  where  a  is  as  de£ined  by  Eq.  4.3.  The  dependent 
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variables  were  v(t)  =  -V(x)  and  c(t)  =  C(x)  +  1/t.  The  starting 
value  for  t  was  normally  taken  as  about  2"^^  3.2x10“^.  The  differ¬ 

ential  equations  were  used  in  the  form  dv/dt  =  M^/(atE)  and  dc/dt 
=  [-1  +  ,  where 

E  =  (l-ct)2  -  ,  (8.3) 

=  v(l-v)  (l-av)t^  -  P4(l-ct)^(v-Pg)  , 

=  (l-ct)2(a  +  ^)  -  t2[(l-v)(l-p22v)  +  p^jV]  , 

and  the  constants  are 

P4  =  aCv+1) 

p  =  2~^ -  (8.4) 

ay (v+1) 

P21  =  (l-a)/Y 

P22  =  a(l  +  1/2v(y-1)) 

P23  =  l/2(Y-l)(l-a). 

The  starting  value  for  v(t)  is  p^,  removing  the  1/t  singularity  in 
dv/dt.  The  starting  value  for  c(t)  is  zero;  it  can  be  determined  by 
substitution  that  the  starting  value  for  M  is  then  just  a,  which 
removes  the  1/t^  singularity  from  dc/dt.  A  little  analysis  shows 
that  v(t)  is  even  and  c(t)  is  odd,  so  there  is  in  fact  no  1/t 
singularity  either. 

The  integration  is  carried  out  until  V  =  -v  and  C  =  c  -  1/t 
match  the  target  functions.  The  interpolated  value  of  x  at  which 
the  match  occurs  is  then  3.  If  it  is  desired  to  tabulate  the  solu¬ 
tion  behind  the  reflected  shock  against  x,  which  runs  from  3  to 
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infinity,  rather  than  against  t,  then  the  parameter  k  can  be  identi- 
fied  as  k  = 

IX.  THE  DENSITY,  AND  THE  MACH  NUMBER  OF  THE  REFLECTED  SHOCK 

Taking  the  initial  density  as  unity,  we  can  use  the  con¬ 

stant  of  the  motion  to  find  that,  in  the  region  between  the  incoming 
shock  and  the  reflected  shock  (call  it  Region  2), 


P2 


=  P2C2C)  = 


Y+1 

'^0 

c 

Y-l 

X 

s' 

(9.1) 


where  =  -1,  say,  C^  =  C(Xq),  and  =  V(Xq),  and  where 


_  l-g 
a (v+1) 


(9.2) 


,  _  2a  (v+l)  _ _ 

b  -  [(v+i)Y  -  (v-i) Ja  - L  • 

For  the  region  behind  the  reflected  shock  (call  it  Region  3) , 
we  can  use  the  jump  condition  to  relate  the  densities  at  3* 


P3(3)  = 


Cy+1)M^ 
(Y-1)m2  +  2 


P2 (3) > 


(9.3) 


where  the  Mach  number,  being  the  magnitude  of  the  ratio  of  fluid 
speed  ahead  of  the  shock,  relative  to  shock  speed,  to  sound  speed 
ahead  of  the  shock,  turns  out  to  be  simply 


M  =  I  (1+V2(3))/C2(3)  !• 


(9.4) 
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Then  we  can  use 


=  p^Cx)  =  p^CP) 


i  C  /  1+V  V 

X  C3(B)  \1+V3(3)  ) 


(9.5) 


REFERENCES 

1.  G.  Guderley,  "Starke  kugelige  und  zylindrische  Verdichtungsstosse 

in  der  Nahe  des  Kugelmittelpunktes  bzw.  der  Zylinderachse , " 
Luftfahrt- Forsch .  302-312  (1942). 

2.  D.  S.  Butler,  "Converging  Spherical  and  Cylindrical  Shocks," 
Armament  Research  Establishment  report  54/54  (1954). 

3.  K.  P.  Stanyukovich,  Unsteady  Motion  of  Continuous  Media 
(Academic  Press,  New  York,  196U) ,  pp ,  bUb-528  . 

4.  Ya,  B.  Zel’dovich  and  Yu.  P.  Raizer,  Physics  of  Shock  Waves  and 
High- temperature  Hydrodynamic  Phenomena  (Academic  Press, 

New  York,  1966)  ,  pp .  785 - 8 0 6 '. 

5.  R.  D.  Richtmyer  and  R.  B.  Lazarus,  "Singularity  Fitting  in 
Hydrodynamical  Calculations,"  Los  Alamos  Scientific  Laboratory 
report  LA-6108-MS  (1975). 


U.S.  GOVERNMENT  PRINTING  OFFICE  1977-777-018/74 


27 


